Analysis of Xenium melanoma spatial data using scanpy and squiqpy¶
Introduction¶
In the rapidly evolving field of genomics, understanding the spatial organization of gene expression within tissues is pivotal for unraveling complex biological processes and disease mechanisms. Traditional single-cell RNA sequencing (scRNA-seq) has provided invaluable insights into cellular diversity, but it often lacks the spatial context necessary to fully comprehend cellular interactions and tissue architecture. Enter spatial transcriptomics—a revolutionary approach that maps gene expression directly within the spatial framework of tissues.
Xenium, developed by 10x Genomics, stands at the forefront of spatial gene expression technologies. Leveraging advanced in situ hybridization techniques, Xenium enables high-resolution, multiplexed profiling of RNA molecules within intact tissue samples. This powerful platform generates rich datasets that capture not only the gene expression profiles of individual cells but also their precise spatial locations and interactions within the tissue microenvironment.
Squipy emerges as a versatile Python library tailored for the analysis of spatial genomics data, including that produced by Xenium. Designed with flexibility and scalability in mind, Squipy offers a comprehensive suite of functionalities for data preprocessing, visualization, spatial clustering, and downstream analyses. Its intuitive interface and compatibility with popular data formats make it an indispensable tool for researchers aiming to decode the spatial dimensions of gene expression.
This tutorial is crafted to guide you through the end-to-end analysis of Xenium data using Squipy. We will begin with data import and preprocessing, delve into quality control and normalization, explore spatial visualization and clustering techniques, and conclude with advanced analyses such as spatially informed differential expression and interaction mapping.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
import squidpy as sq
import spatialdata as sd
from spatialdata_io import xenium
There are a number of curated, high quality spatial datasets that are publicly available on the 10X Genomics website. For this demo, we will download the FFPE Human Skin Primary Dermal Melanoma with 5K Human Pan Tissue and Pathways Panel data. If you have no already downloaded this data, you should go to the "Output and supplemental files" tab and download the "Xenium Output Bundle (full)" file, as well as the supplementary files, Post-Xenium H&E Image (OME-TIFF) and H&E Image Alignment File (CSV).
Below, we load AnnData (h5) file from this Xenium dataset that we have downloaded:
### https://cf.10xgenomics.com/samples/xenium/preview/Xenium_FFPE_Human_Breast_Cancer_Rep1/Xenium_FFPE_Human_Breast_Cancer_Rep1_cell_feature_matrix.h5
#adata = sc.read_10x_h5(filename="../xdd201/Downloads/Xenium_FFPE_Human_Breast_Cancer_Rep1_outs/cell_feature_matrix.h5")
adata = sc.read_10x_h5(filename="/home/thelonius/tmp/spatial_workflows/data/xenium_demo_data/Xeniumranger_V1_hSkin_Melanoma_Add_on_FFPE_outs/cell_feature_matrix.h5")
This displays the size of the matrix of the AnnData object, which can be interpreted as 382 genes in 87499 spatial cells
adata
AnnData object with n_obs × n_vars = 87499 × 382
var: 'gene_ids', 'feature_types', 'genome'
Load the cell index
### https://cf.10xgenomics.com/samples/xenium/preview/Xenium_FFPE_Human_Breast_Cancer_Rep1/Xenium_FFPE_Human_Breast_Cancer_Rep1_cells.csv.gz
#df = pd.read_csv("../xdd201/Downloads/Xenium_FFPE_Human_Breast_Cancer_Rep1_outs/cells.csv")
df = pd.read_csv("/home/thelonius/tmp/spatial_workflows/data/xenium_demo_data/Xeniumranger_V1_hSkin_Melanoma_Add_on_FFPE_outs/cells.csv.gz")
df.set_index(adata.obs_names, inplace=True)
adata.obs = df.copy()
adata.obs
| cell_id | x_centroid | y_centroid | transcript_counts | control_probe_counts | control_codeword_counts | unassigned_codeword_counts | deprecated_codeword_counts | total_counts | cell_area | nucleus_area | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| aaaahdpa-1 | aaaahdpa-1 | 839.648499 | 708.609436 | 319 | 0 | 0 | 0 | 0 | 319 | 222.575164 | 33.009220 |
| aaaakeom-1 | aaaakeom-1 | 836.833252 | 630.847229 | 131 | 0 | 0 | 0 | 0 | 131 | 108.104066 | 23.390938 |
| aaaamacd-1 | aaaamacd-1 | 847.556458 | 634.055603 | 83 | 0 | 0 | 0 | 0 | 83 | 64.392815 | 14.088751 |
| aaabcmaf-1 | aaabcmaf-1 | 828.426575 | 636.969177 | 117 | 0 | 0 | 0 | 0 | 117 | 88.054691 | 6.818594 |
| aaabgold-1 | aaabgold-1 | 837.756653 | 639.850037 | 145 | 0 | 0 | 0 | 0 | 145 | 94.782972 | 21.945938 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| oimimkme-1 | oimimkme-1 | 5896.068848 | 1814.893433 | 217 | 0 | 0 | 0 | 0 | 217 | 79.158909 | 29.712814 |
| oimjcdbg-1 | oimjcdbg-1 | 5897.338379 | 1806.026611 | 184 | 0 | 0 | 0 | 0 | 184 | 61.412502 | 26.235782 |
| oimjcodj-1 | oimjcodj-1 | 5901.827637 | 1746.910278 | 131 | 0 | 0 | 0 | 0 | 131 | 414.805328 | 14.766094 |
| oimjcofe-1 | oimjcofe-1 | 5901.106934 | 1763.095947 | 208 | 0 | 0 | 0 | 0 | 208 | 302.998448 | 20.320313 |
| oimkjkok-1 | oimkjkok-1 | 5895.570801 | 1781.120728 | 104 | 0 | 0 | 1 | 0 | 105 | 136.913755 | 19.191407 |
87499 rows × 11 columns
adata.obsm["spatial"] = adata.obs[["x_centroid", "y_centroid"]].copy().to_numpy()
adata.obsm["spatial"]
array([[ 839.64849854, 708.60943604],
[ 836.83325195, 630.847229 ],
[ 847.55645752, 634.05560303],
...,
[5901.82763672, 1746.91027832],
[5901.10693359, 1763.09594727],
[5895.57080078, 1781.12072754]])
Take a look at the head of the genes in the Xenium panel
print(adata.var['gene_ids'].head())
ACER1 ENSG00000167769 ACHE ENSG00000087085 ACTA2 ENSG00000107796 ADAM12 ENSG00000148848 AHNAK2 ENSG00000185567 Name: gene_ids, dtype: object
gene_names_list = adata.var.index.tolist()
print(gene_names_list)
['ACER1', 'ACHE', 'ACTA2', 'ADAM12', 'AHNAK2', 'AHSP', 'AIF1', 'AIM2', 'AKAP5', 'AKR1C1', 'ALAS2', 'ALOX5AP', 'APCDD1', 'APOD', 'AQP1', 'AQP3', 'AREG', 'ARG1', 'ARHGDIB', 'ASPN', 'ATF3', 'AXL', 'AZGP1', 'BASP1', 'BCAN', 'BCL2A1', 'BHLHE41', 'BIRC3', 'BIRC5', 'C11orf96', 'C15orf48', 'C1QA', 'C1orf54', 'C3', 'C3AR1', 'C5AR1', 'CA1', 'CALCRL', 'CALD1', 'CAPNS2', 'CBFA2T3', 'CCL19', 'CCL2', 'CCL22', 'CCL3', 'CCL4', 'CCL5', 'CCL8', 'CCND1', 'CCND3', 'CCR2', 'CCR7', 'CD14', 'CD163', 'CD1A', 'CD1B', 'CD1C', 'CD207', 'CD27', 'CD274', 'CD276', 'CD34', 'CD3D', 'CD3E', 'CD3G', 'CD4', 'CD40LG', 'CD52', 'CD63', 'CD68', 'CD69', 'CD79A', 'CD83', 'CD8A', 'CD93', 'CDC20', 'CDH2', 'CDH5', 'CDK1', 'CENPF', 'CES1', 'CLDN1', 'CLDN4', 'CLDN5', 'CLEC10A', 'CLEC9A', 'COCH', 'COL17A1', 'COL5A2', 'COL6A1', 'COL6A2', 'COL6A3', 'COMP', 'COTL1', 'CPE', 'CPVL', 'CR2', 'CRABP2', 'CRYAB', 'CSF1R', 'CSPG4', 'CST7', 'CTLA4', 'CTNNB1', 'CTSL', 'CTSW', 'CTSZ', 'CXADR', 'CXCL1', 'CXCL10', 'CXCL11', 'CXCL12', 'CXCL13', 'CXCL2', 'CXCL5', 'CXCR3', 'CXCR4', 'CYB561A3', 'DCT', 'DEFB1', 'DERL3', 'DMKN', 'DSC1', 'DSP', 'DST', 'DUSP2', 'ECSCR', 'EDNRB', 'ENAH', 'ENPP1', 'ENTPD1', 'ETV1', 'FABP4', 'FAM210B', 'FAP', 'FBLN1', 'FBXO7', 'FCER1A', 'FCER1G', 'FCER2', 'FCRL5', 'FGFBP1', 'FKBP8', 'FN1', 'FOXP3', 'FRZB', 'FSCN1', 'G0S2', 'GATA2', 'GEM', 'GIMAP7', 'GNG11', 'GNLY', 'GPR157', 'GPR183', 'GPR4', 'GYPC', 'GZMA', 'GZMB', 'GZMK', 'HAPLN1', 'HAVCR2', 'HBD', 'HBE1', 'HES1', 'HIST1H4C', 'HLA-DPB1', 'HMGB2', 'HMGN2', 'HTRA1', 'ID2', 'ID4', 'IDO1', 'IFI27', 'IFI6', 'IFITM1', 'IFNG', 'IGFBP2', 'IGFBP3', 'IGFBP5', 'IGFBP7', 'IGLL1', 'IL10', 'IL17A', 'IL1B', 'IL2RA', 'IL32', 'IL33', 'IL3RA', 'IL4I1', 'IL6', 'IL7R', 'INHBA', 'INSIG1', 'IRF1', 'IRF4', 'IRF5', 'IRF8', 'ITGA5', 'ITGA6', 'ITGA8', 'ITGAE', 'ITGAX', 'ITGB1', 'ITM2A', 'KDM5B', 'KIT', 'KLRB1', 'KLRD1', 'KLRG1', 'KRT1', 'KRT15', 'KRT16', 'KRT17', 'KRT18', 'KRT2', 'KRT5', 'KRTDAP', 'LAG3', 'LAMB3', 'LCK', 'LENG8', 'LEPR', 'LGALS2', 'LGALS9', 'LINC00636', 'LOR', 'LSAMP', 'LST1', 'LTB', 'LUM', 'LY6D', 'LYPD3', 'LYVE1', 'LYZ', 'MAL2', 'MARCKSL1', 'MEF2C', 'MFAP4', 'MFAP5', 'MGP', 'MIA', 'MIR205HG', 'MITF', 'MKI67', 'MLANA', 'MMP10', 'MMP2', 'MMP27', 'MMP9', 'MMRN1', 'MNDA', 'MPZ', 'MS4A1', 'MT1A', 'MT1F', 'MT1G', 'MT1X', 'MYC', 'MYCT1', 'MYH11', 'MYL9', 'MYOC', 'MZB1', 'NDUFA4L2', 'NFKBIZ', 'NGFR', 'NKG7', 'NLGN4Y', 'NMB', 'NOTCH3', 'NR2F2', 'NR4A2', 'NRP2', 'NSG1', 'NUSAP1', 'PBXIP1', 'PCDH7', 'PDCD1', 'PDCD1LG2', 'PDGFRA', 'PDGFRB', 'PDPN', 'PECAM1', 'PERP', 'PKHD1L1', 'PKIB', 'PLAC9', 'PLP1', 'PLVAP', 'PMEL', 'POSTN', 'PPP1R14A', 'PRF1', 'PROX1', 'PTCH1', 'PTEN', 'PTGDS', 'PTGER3', 'PTGR1', 'PTN', 'PTPRCAP', 'PTTG1', 'QPCT', 'RAMP2', 'RASSF8', 'RGCC', 'RGS1', 'RGS16', 'RGS5', 'RHOV', 'RNASE1', 'RRM2', 'S100A13', 'S100A14', 'S100A4', 'S100B', 'SBSN', 'SCGB1B2P', 'SELE', 'SELL', 'SERPINB2', 'SERPINB5', 'SERPINE2', 'SERPING1', 'SFRP1', 'SFRP2', 'SLC25A37', 'SLC25A39', 'SLC2A3', 'SLC8A1', 'SLPI', 'SNCA', 'SOD3', 'SOSTDC1', 'SOX10', 'SOX15', 'SOX17', 'SPARCL1', 'SPRY1', 'STMN1', 'STRADB', 'SYNPO2', 'TACSTD2', 'TAGLN', 'TBC1D7', 'TCF7', 'TFAP2A', 'TFAP2B', 'TFF3', 'TFPI', 'TFPI2', 'TGFB1', 'THBS2', 'THY1', 'TIGIT', 'TM4SF1', 'TMEM150C', 'TMEM176A', 'TMEM45A', 'TNC', 'TNF', 'TNFRSF17', 'TNFRSF18', 'TOP2A', 'TP63', 'TRPM1', 'TSC22D1', 'TSPAN33', 'TUBB2B', 'TYMS', 'TYR', 'TYROBP', 'TYRP1', 'UBE2C', 'VCAN', 'VEGFA', 'VIM', 'VWF', 'WDFY4', 'WIF1', 'ZNF667-AS1']
It looks like out of all the genes in this Xenium panel, the following genes are likely to be highly expressed in melanoma, but not as highly expressed in healthy tissue. We will visualize some of the following genes later in this workflow: MITF, MLANA, PMEL, TYR, TYRP1, DCT, MIA, and S100B
Calculate quality control metrics¶
Calculate the quality control metrics on the anndata.AnnData using scanpy.pp.calculate_qc_metrics.
### Calculate quality control metrics
sc.pp.calculate_qc_metrics(adata, percent_top=(10, 20, 50, 150), inplace=True)
The percentage of control probes and control codewords can be calculated from adata.obs
### calculate and print the percentage of control probes
cprobes = (adata.obs["control_probe_counts"].sum() / adata.obs["total_counts"].sum() * 100)
print(f"Negative DNA probe count % : {cprobes}")
Negative DNA probe count % : 0.0012551331253257809
### calculate and print the percentage of control codewords
cwords = (adata.obs["control_codeword_counts"].sum() / adata.obs["total_counts"].sum() * 100)
print(f"Negative decoding count % : {cwords}")
Negative decoding count % : 0.003366710030285624
Next we plot the distribution of total transcripts per cell, unique transcripts per cell, area of segmented cells and the ratio of nuclei area to their cells
### plot the distribution of:
fig, axs = plt.subplots(1, 4, figsize=(15, 4))
axs[0].set_title("Total transcripts per cell")
sns.histplot(adata.obs["total_counts"], kde=False, ax=axs[0])
axs[1].set_title("Unique transcripts per cell")
sns.histplot(adata.obs["n_genes_by_counts"], kde=False, ax=axs[1])
axs[2].set_title("Area of segmented cells")
sns.histplot(adata.obs["cell_area"], kde=False, ax=axs[2])
axs[3].set_title("Nucleus ratio")
sns.histplot(adata.obs["nucleus_area"] / adata.obs["cell_area"], kde=False, ax=axs[3])
<Axes: title={'center': 'Nucleus ratio'}, ylabel='Count'>
Filter the cells based on the minimum number of counts required using scanpy.pp.filter_cells. Filter the genes based on the minimum number of cells required with scanpy.pp.filter_genes. The parameters for the both were specified based on the plots above. They were set to filter out the cells and genes with minimum counts and minimum cells respectively.
Other filter criteria might be cell area, DAPI signal or a minimum of unique transcripts.
## Filter the cells and genes
sc.pp.filter_cells(adata, min_counts=10)
sc.pp.filter_genes(adata, min_cells=5)
Observe the size of the adata object after filtering
adata # 167780 × 313
AnnData object with n_obs × n_vars = 87477 × 382
obs: 'cell_id', 'x_centroid', 'y_centroid', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'deprecated_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area', 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'log1p_total_counts', 'pct_counts_in_top_10_genes', 'pct_counts_in_top_20_genes', 'pct_counts_in_top_50_genes', 'pct_counts_in_top_150_genes', 'n_counts'
var: 'gene_ids', 'feature_types', 'genome', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells'
obsm: 'spatial'
Normalize counts per cell using scanpy.pp.normalize_total.
Logarithmize, do principal component analysis, compute a neighborhood graph of the observations using scanpy.pp.log1p, scanpy.pp.pca and scanpy.pp.neighbors respectively.
Use scanpy.tl.umap to embed the neighborhood graph of the data and cluster the cells into subgroups employing scanpy.tl.leiden.
## standard scanpy workflow
sc.pp.normalize_total(adata, inplace=True) # Normalize counts per cell
sc.pp.log1p(adata) # Logarithmize
# sc.pp.highly_variable_genes (adata) # 313 genes
sc.pp.pca(adata) # do principal component analysis
sc.pp.neighbors(adata) # compute a neighborhood graph
sc.tl.umap(adata) # embed the neighborhood graph of the data
sc.tl.leiden(adata, resolution=0.3) # cluster the cells into subgroups
Visualize annotation on UMAP and spatial coordinates¶
Subplot with scatter plot in UMAP (Uniform Manifold Approximation and Projection) basis. The embedded points were colored, respectively, according to the total counts, number of genes by counts, and leiden clusters in each of the subplots. This gives us some idea of what the data looks like.
Visualize cell clusters on UMAP¶
sc.pl.umap(adata, color="leiden", legend_loc='on data', wspace=0.4)
Visualize gene expression on UMAP¶
# sc.pl.umap(adata, color=['ACER1', 'ACTA2'], use_raw=False)
sc.pl.umap(adata, color=['MITF', 'MLANA', 'PMEL', 'TYR', 'TYRP1', 'DCT', 'MIA', 'S100B'], use_raw=False)
Visualize cell clusters on spatial coordinates¶
sq.pl.spatial_scatter(adata, library_id="spatial", shape=None, color=['leiden'], wspace=0.4)
/home/thelonius/anaconda3/envs/spatial_data2/lib/python3.10/site-packages/squidpy/pl/_spatial_utils.py:976: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
Computation of spatial statistics¶
This example shows how to compute centrality scores, given a spatial graph and cell type annotation.
The scores calculated are closeness centrality, degree centrality and clustering coefficient with the following properties:
closeness centrality - measure of how close the group is to other nodes.
clustering coefficient - measure of the degree to which nodes cluster together.
degree centrality - fraction of non-group members connected to group members.
All scores are descriptive statistics of the spatial graph.
This dataset contains Leiden cluster groups’ annotations in adata.obs, which are used for calculation of centrality scores.
Building a spatial neighborhood graph¶
First, we need to compute a connectivity matrix from spatial coordinates to calculate the centrality scores. We can use squidpy.gr.spatial_neighbors for this purpose. We use the coord_type="generic" based on the data and the neighbors are classified with Delaunay triangulation by specifying delaunay=True.
sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)
Compute centrality scores¶
Centrality scores are calculated with squidpy.gr.centrality_scores, with the Leiden groups as clusters.
sq.gr.centrality_scores(adata, cluster_key="leiden")
The results were visualized by plotting the average centrality, closeness centrality, and degree centrality using squidpy.pl.centrality_scores.
sq.pl.centrality_scores(adata, cluster_key="leiden", figsize=(16, 5))
Compute co-occurrence probability¶
This example shows how to compute the co-occurrence probability.
If you want this portion of the workflow to run faster for a demo, it is possible to create a subset of the original Anndata object and work with that (as shown below). You can also work with sdata.tables["subsample"] directly, but to keep verbosity low, we just use the anndata.Anndata object as before.
## how to subsample an AnnData object
'''
adata_subsample = sc.pp.subsample(adata, fraction=0.5, copy=True)
sq.gr.co_occurrence(adata_subsample, cluster_key="leiden")
'''
'\nadata_subsample = sc.pp.subsample(adata, fraction=0.5, copy=True)\nsq.gr.co_occurrence(adata_subsample, cluster_key="leiden")\n'
sq.gr.co_occurrence(adata, cluster_key="leiden")
WARNING: `n_splits` was automatically set to `43` to prevent `87477x87477` distance matrix from being created
0%| | 0/946 [00:00<?, ?/s]
This plot displays the conditional probability of observing one cluster inn the presence of another cluster, where the probability is calculated by observing in the radius size of interest. The score is computed across increasing radii size around each cell in the tissue.
sq.pl.co_occurrence(adata, cluster_key="leiden", clusters="0", figsize=(10, 10))
Neighbors enrichment analysis¶
This example shows how to run the neighbors enrichment analysis routine.
It calculates an enrichment score based on proximity on the connectivity graph of cell clusters. The number of observed events is compared against permutations and a z-score is computed.
This dataset contains cell type annotations in adata.obs which are used for calculation of the neighborhood enrichment. We calculate the neighborhood enrichment score with squidpy.gr.nhood_enrichment.
sq.gr.nhood_enrichment(adata, cluster_key="leiden")
0%| | 0/1000 [00:00<?, ?/s]
And visualize the results with squidpy.pl.nhood_enrichment.
fig, ax = plt.subplots(1, 2, figsize=(13, 7))
sq.pl.nhood_enrichment(adata, cluster_key="leiden", figsize=(8, 8), title="Neighborhood enrichment adata", ax=ax[0])
sq.pl.spatial_scatter(adata, color="leiden", shape=None, size=2, ax=ax[1])
WARNING: Please specify a valid `library_id` or set it permanently in `adata.uns['spatial']`
/home/thelonius/anaconda3/envs/spatial_data2/lib/python3.10/site-packages/squidpy/pl/_spatial_utils.py:976: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored _cax = scatter(
Compute Moran’s I score¶
This example shows how to compute the Moran’s I global spatial auto-correlation statistics.
The Moran’s I global spatial auto-correlation statistics evaluates whether features (i.e. genes) shows a pattern that is clustered, dispersed or random in the tissue are under consideration.
We can compute the Moran’s I score with squidpy.gr.spatial_autocorr and mode = ‘moran’. We first need to compute a spatial graph with squidpy.gr.spatial_neighbors. We will also subset the number of genes to evaluate.
sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)
sq.gr.spatial_autocorr(adata, mode="moran", n_perms=100, n_jobs=1)
adata.uns["moranI"].head(10)
0%| | 0/100 [00:00<?, ?/s]
| I | pval_norm | var_norm | pval_z_sim | pval_sim | var_sim | pval_norm_fdr_bh | pval_z_sim_fdr_bh | pval_sim_fdr_bh | |
|---|---|---|---|---|---|---|---|---|---|
| TYRP1 | 0.907842 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000009 | 0.0 | 0.0 | 0.01014 |
| PMEL | 0.880539 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000008 | 0.0 | 0.0 | 0.01014 |
| KRT16 | 0.800701 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000009 | 0.0 | 0.0 | 0.01014 |
| TYR | 0.800029 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000009 | 0.0 | 0.0 | 0.01014 |
| DSP | 0.784162 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000010 | 0.0 | 0.0 | 0.01014 |
| KRT5 | 0.782448 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000010 | 0.0 | 0.0 | 0.01014 |
| QPCT | 0.781530 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000007 | 0.0 | 0.0 | 0.01014 |
| DMKN | 0.769896 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000007 | 0.0 | 0.0 | 0.01014 |
| KRTDAP | 0.754670 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000009 | 0.0 | 0.0 | 0.01014 |
| KRT17 | 0.753100 | 0.0 | 0.000004 | 0.0 | 0.009901 | 0.000009 | 0.0 | 0.0 | 0.01014 |
We can visualize some of those genes with squidpy.pl.spatial_scatter. We could also pass mode = 'geary' to compute a closely related auto-correlation statistic, Geary’s C. See squidpy.gr.spatial_autocorr for more information.
Let's take a look at some of the genes that have the highest spatial localization:
sq.pl.spatial_scatter(adata, library_id="spatial", shape=None, color=['TYRP1', 'PMEL'], wspace=0.4, size=2, img=False)
sq.pl.spatial_scatter(adata, library_id="spatial", shape=None, color=['TYR', 'QPCT'], wspace=0.4, size=2, img=False)
View histology image¶
Viewing the histology image frequently takes more memory that most local laptops, so it's very important that you have sufficient RAM for the following steps.
# Load the TIFF image into an ImageContainer.
'''
tiff_path = "./data/xenium_demo_data/Xeniumranger_V1_hSkin_Melanoma_Add_on_FFPE_outs/morphology_focus.ome.tiff"
# The ImageContainer supports multi-channel images and allows for efficient I/O.
sq.im.ImageContainer(tiff_path)
# Plot the image.
# The plot() method will render the image using matplotlib.
img.plot()
# Optionally, if you need to further customize or ensure the plot stays visible:
#plt.show()
## For plotting we can also use spatialdata-plot, a package for static plotting of spatialdataobjects.
import spatialdata
#import spatialdata_plot # pip install spatialdata
gene_name = ['ACER1', 'ACTA2']
for name in gene_name:
adata.pl.render_images("morphology_focus").pl.render_shapes(
"cell_circles",
color=name,
table_name="table",
use_raw=False,
).pl.show(
title=f"{name} expression over Morphology image",
coordinate_systems="global",
figsize=(10, 5),
)
'''
'\ntiff_path = "./data/xenium_demo_data/Xeniumranger_V1_hSkin_Melanoma_Add_on_FFPE_outs/morphology_focus.ome.tiff"\n# The ImageContainer supports multi-channel images and allows for efficient I/O.\nsq.im.ImageContainer(tiff_path)\n\n# Plot the image.\n# The plot() method will render the image using matplotlib.\nimg.plot()\n\n# Optionally, if you need to further customize or ensure the plot stays visible:\n#plt.show()\n\n\n\n## For plotting we can also use spatialdata-plot, a package for static plotting of spatialdataobjects.\nimport spatialdata\n#import spatialdata_plot # pip install spatialdata\n\ngene_name = [\'ACER1\', \'ACTA2\']\nfor name in gene_name:\n adata.pl.render_images("morphology_focus").pl.render_shapes(\n "cell_circles",\n color=name,\n table_name="table",\n use_raw=False,\n ).pl.show(\n title=f"{name} expression over Morphology image",\n coordinate_systems="global",\n figsize=(10, 5),\n )\n'
Interactive analysis and visualization with napari-spatialdata¶
Additionally, one can use napari-spatialdata to visualize Xenium data with an interactive GUI. (TODO)
'''
from napari_spatialdata import Interactive
Interactive(sdata)
'''
'\nfrom napari_spatialdata import Interactive\n\nInteractive(sdata)\n'
References¶
https://squidpy.readthedocs.io/en/stable/notebooks/tutorials/tutorial_xenium.html
https://spatialdata.scverse.org/en/latest/tutorials/notebooks/datasets/README.html
https://squidpy.readthedocs.io/en/stable/notebooks/tutorials/tutorial_image_container_zstacks.html
https://spatialdata.scverse.org/projects/napari/en/latest/notebooks/spatialdata.html